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C°  Collocation  Solution  For  Weak-Formulation  Navier-Stokes  Equations 


H.J.  Migliore,  E.T.  Ray,  E.G.  McReynolds 


Portland  State  University 
Portland,  OR  97207 


INTRODUCTION 

Collocation  was  used  to  solve  a  nonlinear  wave  equation  as  applied  to  the 
analysis  of  cable  dynamics  including  hydrodynamic  loading,  surface  excitation, 
and  change  of  length  with  time  (1).  The  collocation  formulation  was  global  in 
nature  and  was  cast  in  terms  of  expressions  for  spatial  derivatives  vrfiich  were 
directly  substituted  in  the  differential  equations  of  motion.  This  approach 
suggested  a  generality  for  application  of  collocation  which  was  not 
particularly  sensitive  to  the  order  of  the  differential  equation,  to 
nonlinearities,  or  to  complex  boundary  equations,  since  the  latter  were 
expressed  explicitly  along  with  the  system  of  residual  equations  at  interior 
collocation  points. 

As  an  extension  to  this  work,  the  modelling  of  real  fluids,  i.e.  solution 
two-dimensional  Navier-Stokes  equations,  were  of  interest  especially  in  terms 
of  interactive  boundary  conditions,  such  as  encountered  in  structure-fluid 
interaction.  In  light  of  addressing  interaction,  a  primative  variable  version 
of  Navier-Stokes  equation  was  adopted  apriori,  rather  than  stream  or  vorticity 
stream  function  approaches  which  might  interfere  with  the  flexibility  required 
along  boundaries  (2).  For  the  sake  of  economy  in  two-D  and  simplicity  (3,4), 
a  penalty  function  formulation  was  employed.  The  solution  of  these 
two-dimensional  Navier-Stokes  equations  was  to  be  based  on  collocation,  using 
the  previous  global  procedure  as  a  starting-off  point.  This  paper  discusses  a 
collocation  solution  for  the  Navier-Stokes  equations  including  the  historical 
evolution  that  led  to  the  developed  procedure:  discrete-Galerkin  collocation 
on  the  interior  and  weak-formulation  Galerkin  using  FEM  shape  functions  on 
element  boundaries. 


GLOBAL  COLLOCATION 

Several  reviews  and  applications  of  collocation  can  be  found  in  Villadsen 
and  Michelsen  (5),  Finlayson  (6,7),  Pinder  (8,9,10),  and  Becker,  Carey  and 
Oden  (11,12),  and  only  an  overview  will  be  given  here.  Collocation's 
theoretical  base  is  Method  of  Weighted  Residual  (MWR)  where  a  trial  solution 
is  substituted  into  the  differential  equation,  resulting  in  a  residual 
expression  reminiscent  of  the  given  differential  equation.  The  residual  is 
minimized  by  weighting  over  the  domain.  In  the  case  of  Galerkin,  the 
weighting  functions  are  the  basis  functions  associated  with  the  trial 
solution,  and  in  the  case  of  collocation,  the  Dirac  delta  has  historically 
been  interpreted  as  the  weighting  function. 


The  optimal  collocation  points  are  not  generally  arbitrary  or  equally 
spaced.  One  approach  (5,6,12)  uses  orthogonal  trial  functions,  such  as 
Legendre,  In  the  Galerkln  Integral.  Gaussian  quadrature  (12)  is  then  employed 
with  two  consequences:  collocation  Is  linked  to  discrete  Galerkln,  and 
collocation  points  are  determined.  By  thus  minimizing  the  residual, 
expressions  are  hypothetically  available  for  determining  coefficients  of  the 
trial  solution,  analogous  to  a  Rayleigh-Ritz  or  Least-Squares  approach.  The 
flexibility  of  collocation,  however,  Is  to  not  directly  consider  R(xi)=G, 
but  rather  to  develop  expressions  for  derivatives  based  on  using  orthogonal 
polynomials  in  a  discrete  Galerkln  approach. 

Starting  with  a  general  trial  solution  with  three  terms. 
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Employing  Lobatto  quadrat ion, 
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Noting  that  *  -1,  x^  ■  0,  x3  *  1;  8  y  »  w2  *  T  ’  w3  8  I 

Taking  the  derivative  of  the  trial  solution 
3  3  P,. 
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Knowing  P|<  and  its  derivative  as  a  function  of  x  and  knowing  a|<  as  a 
function  of  u(xj),  a  derivative  expression  can  be  found  as  a  function  of 
u(xj)  and  x: 

-gj-  *  [A]u  where  an  element  in  matrix  A  is  Aj  *  A { ak ( u ^ ) } 


[2] 


For  a  derivative  at  a  specific  point,  i: 


3u(xi) 
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u(Xj) 


Similarly, 


[A]2u 
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Rather  than  dealing  with  integrals  or  residuals  directly,  Equations  2  and 
3  are  substituted  into  the  PDE,  resulting  in  a  system  of  ODE'S.  Equations  2 
and  3  are  used  in  either  the  x  or  y  direction  for  the  two  dimensional  case. 

For  mixed  derivatives  in  two-D,  Equation  2  is  used  twice,  relative  to  x  and  y. 
In  applying  these  expressions  boundary  conditions  must  be  incorporated  in  the 
form  of  Neumann  or  Dirichlet  for  second  order  problems.  For  a  specified 
derivative,  the  inverse  of  Equation  2  was  used  to  determine  the  value  of 
u(x|<),  X|<  at  the  boundary,  and,  as  in  the  Dirichlet  case,  the  specified 
value  of  u(x|()  could  then  be  used  in  Equations  2  and  3  for  determining 
u(xj)  at  other  Lobatto  points. 

Global  collocation,  using  derivative  expressions  similar  to  Equations  2 
and  3,  was  attempted  on  penalty  function  Navier-Stokes.  The  test  problem  was 
a  driven  cavity  with  conventional  boundary  conditions,  i.e.  no-slip  on  three 
walls.  In  general,  global  collocation  suffered  from  significant  ringing 
whenever  sharp  gradients  were  encountered  and  required  small  time  integration 
steps  vtfien  collocation  points  were  closely  packed.  A  multi-element  approach 
was  needed  in  order  that  one  trial  function  was  not  forced  to  satisfy  sharp 
gradients  of  the  field  variable,  at  boundaries  in  the  case  of  the  driven 
cavity. 


MULTI -ELEMENT  COLLOCATION 

Historically,  elementary  approaches  in  collocation  were  a  direct 
extension  of  global  techniques,  in  that  collocation  points  were  on  Interiors 
of  elements,  and  interelement  points  provide  local  boundary  conditions.  In 
effect,  additional  expressions  were  needed  to  determine  field  variable  values 
at  interelement  points  and  could  be  derived  by  considering  Cl  continuity  (5,6) 
between  elements  for  second  order  problems.  A  slightly  different  approach  was 
successfully  attempted  by  the  authors  on  one-dimensional  advection-diffusion 
where  boundary  expressions  were  developed  from  C°  and  Cl  considerations  (13). 
For  two-D  Navier-Stokes  equations,  a  conventional  C1  collocation  was 
attempted,  and  Interelement  expressions  were  found  to  be  cumbersome,  thus 
Interfering  with  making  changes  In  defining  mesh  geometry  and  complicating  the 
prospects  of  incorporating  distorting  elements. 

Alternatively,  finite  element  functions  have  been  used  as  trial  functions 
in  conjunction  with  collocation.  For  second  order  problems,  C1  Is  required 
between  elements  (12),  and  Hermite  shape  functions  (8,12)  were  used  to 
automatically  satisfy  this  requirement.  Lagrangian  shape  functions  (7,12) 
have  also  been  considered,  but  require  an  additional  set  of  expressions 
related  to  the  C1  requirement. 

In  either  case,  the  C1  requirement  must  be  definitely  enforced,  and  the 
complexity  of  collocation  finite  element  equations  are  Increased  (12)  relative 


to  Galerkin  finite  element  (GFEM).  GFEM  can  be  integrated  by  parts  resulting 
in  less  stringent  C°  between  elements.  However,  the  Dirac  delta  weighting 
function  associated  with  collocation  cannot  be  differentiated  preventing  a 
weak  MWR  formulation.  Because  distorting  elements  would  be  required  for 
future  interaction  work,  simple  element  formulation  was  necessary,  and  a  C° 
approach  was  preferred. 

C°  FEM  COLLOCATION  FOR  NAVIER-STOKES 

Galerkin  condition  applied  to  Navier-Stokes  equations: 

Du 

/„  [  Uf  -  »•!  ]Ndfl  .  0  [«] 

and  N  are  shape  functions  of  the  form  N^=l  at  (x ,y)=(x-f ,y-j )  and  Nj*0 
at  (x,yHxj,yj)*(xi,yi).  Using  7«(NT)  »  NV»T  +  VN*T  , 

Du 

/fl  TE  N  d  a  +  /q  [VN-T  -  v.(NT)  ]do  *  0 


and  using  the  Divergence  Theorem,  the  weak  form  results: 

Du 

/q[^N  +  I-WldO  -  /aa  N  T«n  ds  *  0  [5] 


The  second  integral,  the  "natural"  boundary  condition,  represents  flux  across 
aa  (12).  If  N  are  Lagrangian  finite  element  functions,  C .  will  be  required 
between  elements  and  the  Galerkin  condition  becomes  (12): 


Du 
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a' 

where  e  are  finite  elements  within  3a.  Equation  6  is  the  basis  for  a 
classical  Galerkin-FEM  approach,  and  is  employed  on  element  boundaries  by 
considering  Nj  for  (xi.yj)  only  on  element  boundaries.  Alternatively,  if 
we  selectively  reverse  the  integration-by-parts: 


Du 


l  ♦  i  1^  • I  •  -  /» n  N  I  •  * 0 


[7] 


Equation  7  is  the  basis  for  our  C°  collocation  approach  and  is  used  on  the 
interior  of  elements  by  considering  for  (xi ,y-f )  only  on  Interiors, 
resulting  In: 


Lobatto  quadrature  results  In  the  discrete  Galerkin  expression: 

Ou 

|  U  -  V-D  N|  *  °  [8] 

Du 

R  =  -  v«T  *  0  evaluated  on  interior  [9] 

'  collocation  points 

This  technique  was  implemented  on  Navier-Stokes  equations  and  will  be 
detailed  in  following  sections.  The  specific  form  of  Navier-Stokes  equations 
that  was  solved  will  first  be  reviewed. 


PENALTY  FUNCTION  NAVIER  STOKES 

As  mentioned,  the  problem  of  interest  was  the  solution  to  the  two-D 
Navier-Stokes  equations.  Because  of  compactness  (3)  and  economy  in  two-D  (4), 
a  penalty  function  approach  was  adopted.  The  formulation  will  be  outlined 
below;  further  details  can  be  found  in  Hughes  (3)  and  Oden  (14). 

Navier-Stokes  equations  can  be  written  as  (15,16): 


Du 

p  of  ■  V-<I> 


for  incompressible  fluid,  where 


T  = 


3u 

-P  +  2^ 


L  “1  T5J  +  I?  )  +  2“  J 

For  the  incompressible  case,  continuity  requires  the  additional  expression 
v»us0,  but  for  the  penalty  function  approach,  continuity  is  approximately 
satisfied  by: 

p 

Noting  that  an  expression  for  p  is  available,  the  system  to  be  solved  becomes: 


/  3u  .  3v  1 
“  l-5y+  3*  J 
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[12] 


COLLOCATION  SOLUTION  Of  NAVIER-STOKES  EQUATIONS 

As  discussed,  the  spatial  derivative  expressions  developed  from  our 


5 


collocation  approach  can  be  directly  substituted  in  Equations  10,  11,  and  12. 
The  question  of  boundary  conditions,  whether  local  in  the  sense  of 
interelement  coupling,  or  global  in  the  sense  of  inflow/outflow  requirements 
for  a  physical  problem,  could  not  be  effectively  answered  by  utilizing  Cl 
continuity.  After  considerable  numerical  experimentation,  the  collocation 
approach  that  yielded  reasonable  geometry  flexibility  and  convergence 
properties,  was  collocation  on  the  interior  and  a  finite  element  approach  on 
the  boundary  of  elements  and  of  the  domain. 

A  nine  node  element  was  used  (15,11).  The  central  node,  number  9,  was 
the  only  collocation  point.  Derivative  expressions,  Equations  2  and  3  were 
applied  to  Equations  10  to  12.  Matrix  [A]  has  been  subscripted  and  has  been 
given  the  designation  [B]  when  differentiation  in  the  y  direction  was 
desired.  Repeated  subscript  indicates  summation;  for  Equation  10: 


u9  +  u9A9kuk  +  v9B9kuk  =  ’  p  A9mpm  +  2  p  SAkV  “p  B9mBmkuk  +  p  B9mAmkvk 


and  =  p  Dm£Aikuk  +  p  DmiBlkvk 


m  =  1,2, 


l  *  1,2,  ...  9 


1,2, 


The  derivation  of  the  [o]  matrix  is  detailed  in  Appendix  I. 
Pm,  gives: 

u9  +  u9A9kuk  +  v9B9kuk  =  l?  A9mDmiA«c 


Substituting 


[13] 


+  2p  A9mAtnk  +  p  B9m8mJ  uk  +  £  A9mDm£Blk  +  p  B9nAnJ  vk 

This  time  varying  equation  became  part  of  a  larger  system  that  included 
boundary  equations  relative  to  elements  and  domain  boundaries.  The  discussion 
on  development  of  these  boundary  equations  and  the  time  solution  are  given  in 
the  next  sections. 


FEM  COUPLING  EQUATIONS 

The  Galerkin  FEM  expressions  for  penalty  function  Navier-Stokes  equations 

are: 


/.W#  *  4  *  +  ^  *  ly  l-P  +  2“%)»  Va  ■  * 

/fle  (P  *  )t  ♦?«*  *  0 


[14] 

[15] 

[16] 


fl  denoting  global  minimization  and  flP  denoting  local  minimization  (implying 
C"1  between  elements).  Figure  1  defines  f  and  $  shape  functions  over  an 
element  (11). 


The  utilization  of  Equation  16  is  described  in  Appendix  I.  Equation  14  is 
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a)  Bilinear  used  for  pressure 
FIGURE  1.  Shape  functions 
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b)  Biquadratic  used  for  velocity 


integrated  by  parts  so  as  to  develop  a  weak  formulation  version,  as  was  done 
in  Equation  6: 

3Y  3Y 

y  $  *  »Sf  *  *  V-p  *  *£>  ir +  1$  ♦  £h>e  ■  1  •! 

The  boundary  integral,  8^  ,  will  be  further  discussed  in  Appendix  II. 

Using  Lobatto  quadrature,  not  summing  on  i,  gives: 

wit“i*uiAiAtviBikuk! 


1  +  lE  a1  a  .U.+  -  b1  B  .u.  +  -  b^A  .  v.  j  »  -  B1 
*  p  *  *k  k  p««kk  p*  «k  k J  p  u 


From  Appendix  I,  »  not  summihg  on  i,  then, 

Mi(Ji  *  uiAikuk ♦  'iBikuk>-  -1?  Mlvjk ♦  H  *  1  “j>iB.ki  u» 


-^MXjBjk^W.bXklvk*  i»;  [17] 

where  i  =  1,..8;  «,j,k  s  1 , . . 9;  W-j  are  2-D  weights  defined  in  Appendix  III, 
Similarly,  the  v  velocity  version  of  Equation  17  was  developed  and  these 
equations,  along  with  Equation  13,  and  its  v  velocity  version  were  then  solved 
in  time. 


IMPLEMENTATION 

Assembly 

Equations  13,  17  and  their  "v"  counterparts  were  applied  to  a  master 
element.  The  domain  was  composed  of  master  elements  vrfiich  were  patched 
together.  For  example,  the  equation  governing  a  midside  node,  node  6  on 
element  e  coincident  with  node  8  on  element  e+1  would  be: 


118] 


Du  Du 

L’ m  7VtT->Jdxd*  ^tl/[o  jf  *8  *  Vdlldxdy  *  0 

The  boundary  integral  3  would  only  apply  on  the  domain  boundary  for 
nodes  that  were  part  of  the  system  degrees  of  freedom,  e.g. in-flow,  outflow 
boundaries.  No  slip  boundaries  were  treated  as  constrained  nodes  with 
specified  velocities  and  were  not  part  of  the  degrees  of  freedom.  Looking 
back  at  the  integration  by  parts.  Equation  7,  the  flux  between  elements  does 
not  appear  as  additional  terms  to  Equation  18  because  the  flux  is  generally 
equal  and  opposite  from  each  element.  Only  in  the  case  of  an  applied  flux 
at  an  element  boundary,  J  N  T*n  ds,  would  Equation  17  need  to  be  modified. 

3ne  '  ‘ 


Time  Varying  Equations 

Because  of  the  stiff  terms  associated  with  the  penalty  parameter,  an 
appropriate  ODE  solver  was  required.  LSODE  (17)  was  employed,  a  subroutine 
based  on  Gear's  algorithm,  with  an  Euler  start-up  procedure. 


Appl ications 

The  C°-Collocation  program  was  exercised  on  several  time  varying  flows 
(18).  Results  consist  of  velocity  vector  fields  at  t  =  1.5  seconds,  for  a 
moving  flap  responding  to  fluid  induced  forces  and  spring  restoration  forces. 
Kinematic  viscosity  in  all  cases  was  approximately  that  of  water:  u  =  1.3  x 
10“6  m2/sec.  Figures  2  and  3  show  a  duct  which  is  0.15  m  high  by  0.30  m  long 
composed  of  nine  node  elements  0.03  m  high  and  0.06  m  long.  Figure  4  shows  a 
duct  which  is  0.15  m  high  by  3.5  m  long  composed  of  nine  node  elements  0.03  m 
high  and  length  from  left  to  right  on  the  figure  of  0.95,  .05,  .05,  .05,  and 
2.4  m  respectively.  Each  duct  had  a  0.06  meter  high  pivoting,  rigid  flap  with 
a  specific  density  of  5.0  and  thickness  represented  by  a  line.  A  torsional 
spring  is  located  at  the  base  of  the  flap,  with  spring  stiffness  0.0023 
N/rad.  Boundary  conditions  for  the  velocity  results  shown  in  Figures  2-4  are: 

Top,  bottom,  sides  of  flap:  u  =  0,  v  =  0 
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The  equations  governing  flap  motion  included  fluid  loading  and  in  turn 
provided  changing  geometry  for  the  fluid  boundary.  The  flap  equations  along 
with  the  Navier-Stokes  equations  were  solved  simultaneously.  In  order  to  more 
closely  model  a  structure-fluid  interaction  problem,  elements  were  allowed  to 
distort  as  a  function  of  flap  rotation.  Flap  geometry  was  approximated  by 
using  bilinear  shape  functions,  making  them  subparametric  relative  to  the 
velocity  functions.  This  distortion  required  modification  to  the  derivative 
expressions  as  was  described  for  the  One-D  case.  The  transformation  from 
One-D  weights,  wi ,  to  Two-D  weights,  W-j ,  included  the  effects  of 
distortion. 


FIGURE  2.  Interactive  flap  velocity  field,  inlet  By  =  0 


Figure  4  incorporated  longer  elements  on  the  inlet  and  outlet  than  on  the 
interior,  and  these  variations  seem  to  significantly  affect  stability  and  the 
response  of  the  flap.  Both  Figures  2  and  3  show  instabilities  above  and 
upstream  from  the  top  of  flap,  and  these  instabilities  extend  to  the  inlet. 

In  the  case  of  specifying  tangential  surface  traction,  Bv,  instabilities  are 
manifested  in  inlet  v-component  velocity.  In  the  case  of  specifying  v 
velocity,  instabilities  at  the  inlet  are  manifested  in  u-component  velocity. 
Figure  4  shows  less  instability  due  to  relocating  the  inlet  farther  from  the 
flap. 


0  =  0.04  9  =  0.019  0.10  m/sec  =  - 

FIGURE  3.  Interactive  flap  velocity  field,  inlet  v  *  0 


FIGURE  4.  Interactive  flap  velocity  field,  inlet  By  *  0  and  long  elements 


NODALLY  INTEGRATED  GFEM 

No  additional  attempts  Mere  made  to  control  instabilities  such  as 
incorporating  upwinding  procedures  (19,20),  or  refining  the  mesh  (15),  as  has 
been  done  with  Galerkin  finite  element.  These  modifications  were  not 
made  because  this  particular  C  -col location  approach  performs  the  same  as 
nodally  integrated  GFEM  and  as  a  consequence  will  suffer  from  similar 
instability  problems  and  will  respond  to  similar  fixes.  Regarding 
similarities.  Equation  17  and  its  "v"  counterpart  were  employed  for  all  nine 
nodes  of  the  master  element  with  virtually  identical  results  as  when  using 
Equation  13  for  node  9  and  Equation  17  for  nodes  1  through  8. 


CONCLUSIONS 

Multi-element  collocation  was  successfully  applied  to  the  solution  of  the 
two-dimensional,  unsteady  Navier-Stokes  equations.  In  order  to  take  advantage 
of  the  simplicity  of  C°  between  elements,  a  weak  formulation  version  of  the 
Navier-Stokes  equations  was  desired,  which  precluded  viewing  collocation  as 
based  on  a  Dirac  delta  weighting  function.  The  developed  £r  collocation 
procedure  evaluated  residuals  on  element  interiors  based  on  using  a  quadrated 
version  of  the  MWR  integral.  Element  boundaries  and  coupling  were  treated  as 
in  Galerkin  FEM.  Effectively,  the  developed  C°  collocation  technique  is 
equivalent  to  GFEM  (21).  Because  of  the  normal  (explicit)  form  of  the 
derivative  expressions  and  their  evaluation  at  Lobatto  points,  the  similarity 
between  techniques  is  limited  to  GFEM  that  incorporates  low  order  shape 
functions,  mass  lumping,  and  nodal  integration  (19)  rather  than  Gaussian 
integration.  Because  of  this  equivalency,  no  additional  comparisons  were  made 
between  collocation  and  GFEM  relative  to  finite  difference  or  similar 
techniques  (3,4,15,20).  The  advantages  of  this  C°  collocation  technique  are: 

1.  nodal  values  are  calculated  without  using  additional  expressions  for 
quadrature  and  interpolation 


2.  the  system  of  time  equations  are  in  normal  (explicit)  form 
simplifying  the  incorporation  of  other  governing  equations  as 
encountered  in  an  interaction  problem 

3.  shape  functions  are  not  differentiated;  rather  derivative  expressions 
based  on  Lobatto  quadrature  and  Legendre  polynomials  were  developed 
facilitating  the  mixing  of  shape  functions  as  is  required  with 
penalty  function  approach,  and  facilitating  the  quadration  of 
integral  equations:  both  interior  GFEM  Navier-Stokes  equations  and 
boundary  integrals 

4.  derivative  expressions  are  used  directly  in  residual  equations, 
specified  derivatives  and  other  derivative  equations. 

These  simplifying  characteristics  proved  to  be  significant  in  attempting 
to  model  and  solve  fluid  interaction. 
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APPENDIX  I.  Theory  for  MD“  Matrix 

Galerkin  integral  for  penalty  expression  (12,14): 

/  eiP  +  x  Ggjf  +  -gp)}  ♦id«e  *  0  1  *  1....4,  bilinear  shape  functions 

Using  Lobatto  quadrature  gives: 

WjPj4ji  +  XWj^Ajkuk  +  Bjkvk^ji  *  0  J  ■  Lobatto  points 

where  &e  domain  is  restricted  l  subdomain  element,  e. 


pj  =  *jkpk  j 

wjVk*j1  ■  -xWk  *  Vk)fj) 

M,kpk =  -1VjiVk  -  AWj*ji8okvk  "here  M’k=  Vu+jk 
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l.e.  ♦  evaluated  at  point  j 
k 


Using  repeated  subscripts  to  indicate  summation: 


"-Ak 

Where  m 


■  -^siVjiVk  -  J"$iwj*jiBjkvk 
01.Mik  =  ,  and  ig  is  Kronecker  delta 
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P«  3  *«gPg  3  ’x*«0n,0iWj*jiAjkuk  '  x*-0m0l Wj  i  Bjkvk 

D«j  3  ♦«0m0iwj*ji 


APPENDIX  II.  Boundary  Conditions 

Domain  boundary  conditions  can  be  a  combination  of  given  velocities  or 
specified  surface  traction  force.  Velocities  are  defined  on  that  portion  of 
the  boundary  3%,  and  surface  tractions  are  defined  on  3%,  such  that: 

3Q  =  3%  U  3Qu  3%  n  3QU  =  0 

Looking  at  the  boundary  integral,  3;  n  is  outward  normal  to  3%: 

6u s  ^30,  MTn  i,n. +  Ti2  ral ds 

for  n  oriented  along  the  x-axis, 

•1  -  /30hNi  T“ ds 

specifying  the  normal  traction  component  for  the  inflow  boundary,  Tu*k: 

6i  ■  k  ds  ■  «i k 

For  the  outflow  boundary,  the  normal  component  will  be  specified  as  zero, 

Tii  «  0: 


Considering  the  boundary  integral  for  the  second  Navier-Stokes  equation: 

*  /aa  +  T^il-nds 

Again  with  n  oriented  along  x  axis. 


.1 


v  *  ^3«hNi  Tzi  ds  *  V2* 


when  the  shear  traction  component  is  specified  as  zero,  T21  *  0: 


APPENDIX  III.  Integration  Weights 


The  master  element  is  defined  in  terms  of  coordinates,  £  and  n,  on  a 
global  Cartesian  coordinate  framework,  x,y.  The  integration  weights  must  also 
be  modified  as  follows: 


/a  f(x,y)dxdy  *  _i F(  E,n)  |  j|d£dn 
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where 


3x  3y 


3x  3y 
’  "5n  tfc 


where  w1.w2.w3  are  one-dimensional  Lobatto  weights.  Thus,  integration  of  any 
function  f(x,y)  e.g.  u,v,  has  an  associated  Lobatto  quadrature  expression  in 
terms  of  modified  weights  W1...W9  and  functions  F(E,n)  on  the  master  element. 
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NOMENCLATURE 


A,  A-j j  x-derivative  matrix 

B, Bjj  y-derivative  matrix 

a|<  coefficients  in  general  solution 
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equal  to  ,  accounts  for 


3*t 

aa^ ,  accounts  ror  — gy 

zero'th  derivative  continuity 
first  derivative  continuity 
total  derivative 
elements 

x,y  unit  vectors 

subscripts  designating  variables 

determinant  of  Jacobian  matrix 


evaluated  at  collocation  points 


wj*ij*jk 
5*  ;  m  =  M-1 

normal  vector  to  3ft 
shape  function 
Legendre  polynomials 
pressure 
residual 

integration  variable  associated  with  3ft 
surface  traction  matrix 
components  of  T 
time 

one  dimensional  field  variable,  such  as  velocity 
velocity  vector 

velocity  components  in  x,y  direction 
Lobatto  quadrature  weights 
master  element  integration  weights 
global,  rectangular  coordinates 
Lobatto  quadrature  points 
boundary  integral 

Kronecker  delta 

angle  of  flap  rotation,  measure  clockwise  from  vertical,  in  radians 


0  angular  velocity  of  flap,  in  radians/second 

X  large  number,  penalty  parameter 

u  viscosity 

u  kinematic  viscosity 

p  fluid  density 

5,n  element  coordinates 

bilinear  shape  functions 
biquadratic  shape  functions 
ft  spatial  domain 

3ft  boundary  of  domain,  ft 

v  Two-D  grad  operator 

v«T  vector  surface  force 
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